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Introduction 


The Space Exploration Initiative is considering use of low thrust 
(nuclear electric, solar electric) and intermediate thrust (nuclear thermal) 
propulsion systems for transfer to Mars and back. Due to the duration of 
such a mission, a low thrust minimum-fuel solution is of interest; a savings 
of fuel can be substantial if one allows the propulsion system to be turned 
off and back on. This switching of the propulsion system helps distinguish 
the minimum-fuel problem from the well-known minimum-time problem. 

Optimal orbit transfers are also of interest to the development of a 
guidance system for orbital maneuvering vehicles which will be needed, for 
example, to deliver cargoes to Space Station Freedom. 

The problem of optimizing trajectories for an orbit-to-orbit transfer 
with minimum-fuel expenditure using a low thrust propulsion system was 
discussed in (1), but the code, SECKSPOT (2), which was used is incapable 
of handling a general minimum propellant problem. 


Analysis 


To avoid the singularities that may occur when using the classical 
elements a, e, i, <b, Q, (semi-major axis, eccentricity, angle of inclination, 
argument of perigee, longitude of the ascending node) the equinoctial elements 
a, h, k, p, and q are used to define the spacecraft's state. The sixth state 
element is the eccentric longitude F, which represents the angular position 
of the satellite in its orbit. Thus, the spacecraft's state, z, is 

T 

z ■ (a, h, k, p, q) , or in classical elements, 

z ■ (a, esin(o) + Q) , ecos (o> + £1), tan (i/2) sin (Q) , tan (i/2) cos (Q) ) T , 


and 


[ 1 ] 


F - E + e» + O. 


The convention of denoting vectors in bold print 
the circumflex a is used. 


and unit vectors with 


It follows that the equations of motion in terms 
elements are 


of these equinoctial 


. A 

z » 2P Mu / mc f 


[ 2 ] 


m - -2P / c 2 , 

where P (k 0) " power due to thrusters, m (> 0) ■ mass of the spacecraft, 
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c (> 0) * jet exhaust speed, M » 5 x 3 matrix representing the partial 
derivative of z with respect to the velocity vector r, and the unit vector 

A 

in the direction of thrust is u, which is also the control for this optimi- 
zation problem. 

A 

Applying the maximum principle to the optimal pair (u, z ) , it is 
clear that one should maximize the Hamiltonian function, H, with respect 
to the thrust direction at each point along the optimal trajectory of 
the transfer to obtain the minimum-fuel solution^ The Hamiltonian function 
is defined in terms of the costate (adjoint) variables and X m so that 


• * 
z = aH / 91^, m ■ 3H / dX m , 

X ■ - dH / ds, X m - - 9H / d\ m . 


[3] 


where 

H<x, m, X x , X , u) • ZB / me M u - 2P / c 2 X ffl . [ 

Hence, the Hamiltonian H in [4] is maximized by aligning u in the direction 

of the primer vector, \ ^ M. Let \ denote the magnitude of this primer 

vector. Then define the switching function, o, by o - X - mX m / c. With 
this definition, H can be rewritten as H ■ (2P / mc)o. 


At any point along the trajectory for which a, and hence H, is nega- 
tive, H can be maximized by talcing P - 0, which amounts to turning the 
propulsion system off; while if c, and hence H, is positive, then allow the 
propulsion system to be used at full power in the optimal direction of thrust 

u. So o < 0 determines a "coast phase", while o > 0 yields a "thrust phase". 
The points at which a • 0 are called switch points. 

Around each orbit the equinoctial elements vary "slowly" and can be 
held constant; thus only the "fast" variable F remains. However, in this 
situation the actual position in an orbit of the spacecraft is not of interest, 
so the dependence on this variable is removed by integration. Using Kepler's 
equation which relates the time t to the mean anomaly and hence the eccentric 
anomaly, one cam perform averaging over an orbit of some specified period by 
changing the variable of integration from time t to the eccentric longitude F. 
Thus the Hamiltonian and hence the equations of motion can be averaged; note 
that since the Hamiltonian is zero for coast phases, one needs only integrate 
over thrust phases. In order to perform this integration using a Gaussian 
quadrature method, it is necessary to predetermine the switch points along 
the trajectory and integrate only over thrust periods, which will begin and 
end with a switch point. These averaged derivatives are then solved using 
a Runge-Kutta method. 
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To find th© solution it is necessary to determine the initial values 
for the costate variables X such that the boundary conditions are satisfied; 
that is, the desired final ftate is reached. The transversality conditions 
at the final time (which is unspecified) t f , are * a, where o is a 

vector of parameters and X m (t f ) ■ 1. 


Numerical Algorithm 

The minimum-fuel solution is found by an iterative method which 
determines the initial costate values Xz(t Q ) and X m (t Q ) such the boundary 

conditions a - a f ■ 0, h - h f * 0, k - k f - 0, p - p. ■ 0, q - q f » 0 and 
X m j - 1 ■ 0 (and f - f f if not averaging) are satisfied. This is a 
system of nonlinear equations in terms of whose solution can be 

found by a secant method. For an initial guess, an optimal trajectory 
is found in the following manner: call the Runge-Kutta integrator, which 

calls RKFCT which calculates the derivatives . RKFCT calls SWCHPTS to find 
the thrust periods over and orbit, by finding the switchpoits and calls INTEG to 
prepare the integrands for quadrature if AVGTST determines that the variables 
are changing slowly enough, and if not averaging the partial derivatives 
of the equations of motion are calculated; INTEG calls QUAD 4 for integration. 
QUAD4 calls FCT, which calculates the Hamiltonian, and calls EVALMP and 
PRIMER to find the thrust direction. From this Runge-Kutta integrator, 
a new state and costate is found. The stopping criterion used is the size of 
the semi -major axis. 


Conclusions 

It appears that one of the difficulties involved with this 
problem is convergence of the method considering the sensitivity to 
slight perturbations in the initial costate values. Another problem 
is the accurate determination of the switch points for the case of 
averaging and also not averaging. 
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